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Quasiparticle excitations can compromise the 
performance of superconducting devices, caus¬ 
ing high frequency dissipation, decoherence in 
Josephson qubits GHe], and braiding errors in 
proposed Majorana-based topological quantum 
computers GH2J- Quasiparticle dynamics have 
been studied in detail in metallic superconduc¬ 
tors jToffTH but remain relatively unexplored in 
semiconductor-superconductor structures, which 
are now being intensely pursued in the context 
of topological superconductivity. To this end, 
we introduce a new physical system comprised of 
a gate-confined semiconductor nanowire with an 
epitaxially grown superconductor layer, yielding 
an isolated, proximitized nanowire segment. We 
identify Andreev-like bound states in the semi¬ 
conductor via bias spectroscopy, determine the 
characteristic temperatures and magnetic fields 
for quasiparticle excitations, and extract a par¬ 
ity lifetime (poisoning time) of the bound state 
in the semiconductor exceeding 10 ms. 

Semiconductor-superconductor hybrids have been in¬ 
vestigated for many years [TSinU] . but recently have re¬ 
ceived renewed interest in the context of topological su¬ 
perconductivity, motivated by the realization that com¬ 
bining spin-orbit interaction, Zeeman splitting and prox¬ 
imity coupling to a conventional s-wave superconductor 
provides the necessary ingredients to create Majorana 
modes at the ends of a one-dimensional (ID) wire. Such 
modes are expected to show nonabelian statistics, allow¬ 
ing, in principle, topological encoding of quantum infor¬ 
mation [[20H221 among other interesting effects [231 EH ■ 

Transport experiments on semiconductor nanowires 
proximitized by a grounded superconductor have recently 
revealed characteristic features of Majorana modes UB- 
128] . Semiconductor quantum dots with superconduct¬ 
ing leads have also been explored experimentally [25H52] , 
and have been proposed as a basis for Majorana chains 
[33fl35] . Here, we expand the geometries investigated 
in this context by creating an isolated semiconductor- 
supercondutor hybrid quantum dot (HQD) connected to 
normal leads. The device forms the basis of an isolated 
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Majorana system with protected total parity, where both 
the semiconductor nanowire and the metallic supercon¬ 
ductor are mesoscopic EMEU- 

The measured device consists of an InAs nanowire with 
epitaxial superconducting A1 on two facets of the hexag¬ 
onal wire, with Au ohmic contacts (Figs. la,b). Four de¬ 
vices showing similar behavior have been measured. The 
InAs nanowire was grown without stacking faults using 
molecular beam epitaxy with A1 deposited in situ to en¬ 
sure high-quality proximity effect 03 S3- Differential 
conductance, g , was measured in a dilution refrigerator 
with base electron temperature T ~ 50 mK using stan¬ 
dard ac lock-in techniques. Local side gates, patterned 
with electron beam lithography, and a global back gate 
were adjusted to form an Al-InAs HQD in the Coulomb 
blockade regime, with gate-controlled weak tunneling to 
the leads. The lower right gate, Hr, was used to tune the 
occupation of the dot, with a linear compensation from 
the lower left gate, Vl, to keep tunneling to the leads 
symmetric. We parameterize this with a single effective 

a 


FIG. 1: Nanowire-based hybrid quantum dot. a, Scanning 
electron micrograph of the reported device, consisting of an InAs 
nanowire (gray) with segment of epitaxial A1 on two facets (blue) 
and Ti/Au contacts and side gates (yellow) on a doped silicon sub¬ 
strate. b, Device schematic and measurement setup, showing ori¬ 
entation of magnetic field, B. c, Differential conductance, g , as 
a function of effective gate voltage, Vq, and source-drain voltage, 
Vsd, at B = 0. Even (e) and odd (o) occupied Coulomb valleys 
labeled. 
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gate voltage, Vq (see Supplement). 

Differential conductance as a function of Vg and 
source-drain bias, Vsd, reveals a series of Coulomb dia¬ 
monds, corresponding to incremental single-charge states 
of the HQD (Fig. lc). While conductance features at 
high bias are essentially identical in each diamond, at 
low bias, Vsd < 0.2 mV, a distinctive even-odd pattern 
of left- and right-facing conductance features is observed. 
This results in an even-odd alternation of Coulomb block¬ 
ade peak spacings at zero bias, similar to even-odd spac- 
ings seen in metallic superconductors mum. However, 
the parity-dependent reversing pattern of subgap fea¬ 
tures at nonzero bias has not been reported before, to 
our knowledge. The even-odd pattern indicates that a 
parity-sensitive bound state is being filled and emptied 
as electrons are added to the HQD. 

Measured charging energy, Ec = 1.1 meV, and su¬ 
perconducting gap, A = 180 /J.eV, satisfy the condition 
(A < Ec) for single electron charging [42l 143] . Differ¬ 
ential conductance at low bias occurs in a series of nar¬ 
row features symmetric about zero bias, suggesting trans¬ 
port through an Andreev-like bound state, with negative 
differential conductance (NDC) observed at the border 
of odd diamonds. NDC arises from slow quasiparticle 
escape, as discussed below, similar to current-blocking 
seen in metallic superconducting islands in the opposite 
regime, A > Ec (HI 115] . 

To gain quantitative understanding of these features, 
we model transport through a single Andreev bound state 
in the InAs plus a Bardeen-Cooper-Schriffer (BCS) con¬ 
tinuum in the Al. The model assumes symmetric cou¬ 
pling of both the bound state and continuum to the 
leads, motivated by the observed symmetry in Vsd of 
the Coulomb diamonds. Transition rates were calculated 
from Fermi’s golden rule and a steady-state Pauli mas¬ 
ter equation was solved for state occupancies. Conduc¬ 
tance was then calculated from occupancies and transi¬ 
tion rates (see Supplement). 

Measured and model conductances are compared in 
Figs. 2a,b. The coupling of the bound state to each 
lead, noting the near-symmetry of the diamonds, was es¬ 
timated to be r 0 = 0.5 GHz, based on zero-bias conduc¬ 
tance (Fig. 2d). The energy of the discrete state, E 0 = 
58 /zeV at zero magnetic field, was measured using finite 
bias spectroscopy (Fig. 2e). The normal-state conduc¬ 
tance from each lead to the continuum, g&\ = 0.15 e 2 //i, 
was estimated by comparing Coulomb blockaded trans¬ 
port features in the high bias regime (Vsd = 0.4 mV). 
The superconducting gap, A = 180 /xeV, was found 
from the onset of NDC, which is expected to occur at 
eVsD = A — E 0 (Fig. 2f). While the rate model shows 
good agreement with experimental data, some features 
are not captured, including broadening at high bias, with 
greater broadening correlated with weaker NDC, and 
peak-to-peak fluctuations in the slope of the NDC fea¬ 
ture. These features may be related to heating or cotun- 
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FIG. 2: Subgap bias spectroscopy, experiment and model. 

a, Experimental differential conductance, g, as a function of gate 
voltage Vg and source-drain Vsd, shows characteristic pattern 
including negative differential conductivity (NDC). b, Transport 
model of a. vg = oV g up to an offset, where a is the gate lever 
arm. Axis units are A/e = 180 pV, where A is the supercon¬ 
ducting gap. See text for model parameters, c, Source and drain 
(gold) chemical potentials align with the middle of the gap in the 
HQD density of states. No transport occurs due to the presence 
of superconductivity, d, Discrete state in resonance with the leads 
at zero bias. Transport occurs through single quasiparticle states, 
e, Discrete state in resonance with the leads at high bias. Trans¬ 
port occurs through single and double (particle-hole) quasiparticle 
states, f, Discrete state and BCS continuum in the bias window. 
Transport is blocked when a quasiparticle is in the continuum, re¬ 
sulting in NDC. 


neling, not accounted for by the model. 

The observation of negative differential conductance 
places a bound on the relaxation rate of a single quasi¬ 
particle in the HQD from the continuum (in the Al) to 
the bound state (in the InAs nanowire). Negative differ¬ 
ential conductance arises when an electron tunnels into 
the weakly coupled BCS continuum, blockading trans¬ 
port until it exits via the lead. The blocking condition is 
shown for a hole-like excitation in Fig. 2f. Unblocking oc¬ 
curs when the quasiparticle relaxes into the bound state, 
followed by a fast escape to the leads. NDC thus indi¬ 
cates a long quasiparticle relaxation time, r qp , from the 
continuum to the bound state. Using independently de¬ 
termined parameters, the observed NDC is only compat¬ 
ible with the model when r qp > 0.1 /is (see Supplement). 
This bound on r qp is used below to similarly constrain 
the characteristic poisoning time for the bound state. 

Turning our attention to the even-odd structure of 
zero-bias Coulomb peaks (Figs. 3a,b), we observed con¬ 
sistent large-small peak spacings (Fig. 3), associating the 
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eludes thermal quasiparticle excitations follows earlier 
treatments [40] El 33] > including a discrete subgap state 
as well as the BCS continuum [44] (Fig. 3d). Even-odd 
°<q 5 peak spacing difference, S e — S Q , depends on the differ- 
ence of free energies, 

Se - So = — (Fo - F e ) , (1) 

cue 

where a is the (dimensionless) gate lever arm. The free 
energy difference, written in terms of the ratio of parti¬ 
tion functions, 

F 0 -F e = -k B T\n(^j , (2) 

depends on D(E), the density of states of the HQD, 


FIG. 3: Even-odd Coulomb peak spacings a, Measured zero- 
bias conductance, g, versus gate voltage, Vq , at temperature T ~ 
50 mK, and magnetic field B = 0. b, Peak spacing, S, versus gate 
voltage. Black points show spacings from a calculated using the 
peak centroid (first moment), red points T = 350 mK and B = 0, 
purple points B = 150 mT and T ~ 50 mK. c, Right-most peaks in 
a. Peak maxima (A) and centroids (■) are marked, d, Free energy, 
F, at T = 0 versus gate-induced charge, N , for different HQD 
occupations, where N = CVq/c up to an offset and C is the gate 
capacitance. Parabola intersection points are indicated by circles, 
corresponding to Coulomb peaks. BCS continuum (shaded), shown 
for odd occupancy. Odd Coulomb diamonds carry an energy offset 
Eq for quasiparticle occupation of the sub gap state, resulting in a 
difference in spacing for even and odd diamonds. 


larger spacings with even occupation, as expected the¬ 
oretically [12] 32) and already evident in Fig. 1. Oc¬ 
casional even-odd parity reversals on the timescale of 
hours were observed in some devices, similar to what 
is seen in metallic devices [14]. Peak spacing alterna¬ 
tion disappears at higher magnetic fields, B , consistent 
with the superconducting-to-normal transition, and also 
disappears at elevated temperature, T > 0.4 K, signif¬ 
icantly below the superconducting critical temperature, 
T c ~ 1 K. The temperature dependence is consistent 
with similar behavior seen in metallic structures mm, 
and can be understood as the result of thermal activation 
of quasiparticles within the HQD with fixed total charge. 

As seen in Fig. 3c, individual Coulomb peaks are asym¬ 
metric in shape, with their centroids (first moments) on 
the even sides of the peak maxima. Note that the asym¬ 
metry leads to higher near-peak conductance in even 
valleys, the opposite of the Kondo effect. The asym¬ 
metric shape is most pronounced at low temperature, 
T < 0.15 K, and decreases with increasing magnetic field. 
The degree of asymmetry is not predicted by the rate 
model, even taking into account the known small asym¬ 
metry due to spin degeneracy [461 . In the analysis below, 
we consider peak positions defined both by peak maxima 
and centroids. 

A model of even-odd Coulomb peak spacing that in- 


dE D(E) lncoth[E/(2fc B T)], (3) 

where D(E) consists of one subgap state and the contin¬ 
uum. For A k B T, this can be written 

F 0 - F e « -k B Tln{N eS e- A ' kBT + 2e~ E °' kBT ), (4) 

where N e g = pAiV\/2nk B T A is the effective number of 
continuum states for A1 volume, V , and normal density 
of states pa\ [ID] |4T] (see Supplement). 

Within this model, one can identify a characteristic 
temperature, T* ~ A/[fes ln(IVeff)], less than the gap, 
above which even-odd peak spacing alternation is ex¬ 
pected to disappear. Note in this expression N e g it¬ 
self depends on T, and also that T* does not depend 
on the bound state energy, E 0 . A second (lower) char¬ 
acteristic temperature, T** ~ (A — E 0 )/[k B ln(IV e ff/2)], 
which does depend on Eq, is where the even-odd alter¬ 
nation is affected by the bound state, leading to satu¬ 
ration at low temperature [40] [44] . For a spin-resolved 
zero-energy (E 0 = 0) bound state—the case for unsplit 
Majorana zero modes—these characteristic temperatures 
coincide and even-odd structure vanishes, as pointed out 
in Ref. [3Bj. In the opposite case, where the bound state 
reaches the continuum ( Eq = A), the saturation temper¬ 
ature vanishes, T** = 0, and the metallic result with no 
bound state is recovered mm- 

Experimentally, the average even-odd peak spacing dif¬ 
ference, (S e — S 0 ), was determined by averaging over a set 
of 24 consecutive Coulomb peak spacings, including those 
shown in Fig. 3, at each temperature. Figure 4 shows 
even-odd peak spacing difference appearing abruptly at 
T onset ~ 0.4 K, and saturating at T sat ~ 0.2 K, with a 
saturation amplitude near the value expected from the 
measured bound state energy, 4 Vo = 4E 0 /(ae). Fig¬ 
ure 4 shows good agreement between experiment and the 
model, Eq. ([I]), using a density of states determined inde¬ 
pendently from data in Fig. 2, with V = 7.4 x 10 4 nrn 3 as 
a fit parameter, consistent with the micrograph (Fig. la), 
and pai = 23 eV -1 nm -3 [T4] . 
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FIG. 4: Temperature and magnetic field dependence of the even-odd peak spacings. Average even-odd spacing difference, 
(S e — S 0 ), versus temperature, T. Spacing between peak maxima (triangle) and centroids (square) are shown. Spacing expected from 
lower Zeeman-split bound state, 4Vo(-B) = 4:Eo(B)/(ae), indicated on right axis. Quasiparticle activation temperature, T*, and crossover 
temperature, T** , indicated on top axis. Solid curve is Eq. 0 with a HQD density of states measured from Fig. 2 (A = 180 //eV, 
Eg = 58 //eV. a = 0.013), and the fitted aluminum volume, V = 7.4 X 10 4 nm 3 . Dotted curve includes a discrete state broadening, 
7 = 50 neV, fit to the centroid data. Left inset: Same as main, but at B = 40, 80,100,150 mT, from top to bottom. Curves are fit to two 
shared parameters: g-factor, g = 6, and superconducting critical field, B c = 120 mT, with other parameters fixed from main figure. Right 
inset: Representative Coulomb peaks showing even (S e ) and odd ( S a ) spacings. 


The asymmetric peak shape complicates measurement 
of even-odd spacings, as one can either use the centroids 
or maxima to measure spacings, the two methods giving 
different results. Larger peak tails on the even valley side 
cause the centroids to be more regularly spaced than the 
maxima. This is evident in Fig. 4, where the centroid 
method shows a decreasing peak spacing difference at 
low temperature, while with the maximum method the 
spacing remains flat. The thermal model of S e — S Q can 
also show a decrease at low temperature if broadening of 
the bound state is included (See Methods). We do not 
understand at present if the low temperature decrease in 
the centroid data is related to the decrease seen in the 
model when broadening is included. It is worth noting, 
however, that the fit to the centroid data gives a broaden¬ 
ing 7 = 50 neV, reasonably close to the value estimated 
from the lead couplings, (hT 0 ) 2 /A = 20 neV. 

Applied magnetic field (direction shown in Fig. lb) 
reduces the characteristic temperatures T onset , T sat , and 
saturation amplitudes. Field dependence is modeled by 
including Zeeman splitting of the bound state and orbital 
reduction of the gap and bound state energy, taking the 
g-factor and critical magnetic field as two fit parameters 
applied to all data sets. The fit value g = 6 lies within the 


typical range for InAs nanowires 03 SSI, supporting our 
interpretation of the bound state residing in the InAs. 
The fit value of critical field, B c = 120 mT, is typical for 
this geometry. 

Good agreement between the peak spacing data and 
the thermodynamic model (Fig. 4) suggests that the 
number of thermally activated quasiparticles obeys equi¬ 
librium statistics, N eq (T ) = N 2 s e~ 2A ^ kBT (see Supple¬ 
ment for derivation). Saturation caused by the bound 
state means that even-odd amplitude loses sensitivity as 
a quasiparticle detector below T sat . We therefore take 
N eq (T sat ) ~ 10 -5 (for T sat ~ 0.2 K) as an upper bound 
for the number of quasiparticles at temperatures below 
T sat . The corresponding upper bound of the quasiparticle 
fraction, x qp = N eq (T sat ) / (p A \V A) ~ 10~ 8 , is compara¬ 
ble to values in the recent literature, 10 -5 — 10 -8 , for 
metallic superconducting junctions and qubits j3l-(6} H~3l. 

We now discuss the implications of our measurements 
for determining the poisoning time, r p , of the bound 
state. For the present geometry, the dominant source 
of poisoning of the bound state is not tunneling of elec¬ 
trons from the leads, which is negligible in the strongly 
blockaded regime, but is rather the continuum in the 
strongly-coupled Al, within the isolated structure itself. 
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Theoretical estimates |Hj suggest an inverse relationship 
between t p and the number of available quasiparticles, 
with a proportionality that depends on system details. 
Taking the bound on single quasiparticle relaxation time 
from the continuum into the bound state, r qp > 0.1 /zs, 
from above, as the poisoning time when a single quasipar¬ 
ticle is present, we estimate t p by scaling for the actual 
number of quasiparticles in equilibrium, N eq , giving a 
poisoning time r p = T qp /N eq >10 ms. 

We expect t p to depend weakly on the bound state 
energy for low-energy bound states mm [50] . includ¬ 
ing for Majorana zero modes at Eq = 0. Device geom¬ 
etry may somewhat alter the number of quasiparticles 
available to relax into the bound state, i.e. by chang¬ 
ing iV e ff, but any increase can be compensated by ex¬ 
ponentially small decreases in the quasiparticle temper¬ 
ature. The long poisoning time obtained here suggests 
that a large number of braiding operations in Majorana 
systems should be readily achievable within the relevant 
time scale. 

Methods 

Sample preparation: InAs nanowires were grown in 
the [001] direction with wurzite crystal structure with A1 
epitaxially matched to [111] on two of the six {1100} side- 
facets. They were then deposited randomly onto a doped 
silicon substrate with 100 nm of thermal oxide. Electron- 
beam lithographically patterned wet etch of the epitaxial 
A1 shell (Transene A1 Etchant D, 55 C, 10 s) resulted in a 
submicron A1 segment (310 nm, Fig. 00- Ti/Au (5/100 
nm) ohmic contacts were deposited on the ends following 
in situ At milling (1 mTorr, 300 V, 75 s), with side gates 
deposited in the same step. For the present device, the 
end of the upper left gate broke off during processing. 
However, the device could be tuned well without it. 

Master equations: The master equations (used for 
Fig. lb) consider states with fixed total parity, composed 
of the combined parity of quasiparticles in the thermal- 
ized continuum and the 0, 1, or 2 quasiparticles in the 
bound state (see Supplement). 

Free energy model: Even and odd partition functions 
in Eq. 2, F 0 — F e = —kBThi(Z 0 /Z e ), can be written as 
sums of Boltzmann factors over respectively odd and even 
occupancies of the isolated island. For even-occupancy, 

Z e = l + ^2e- Ei/kBT e - Ej / kB T +..., (5) 

where the first term stands for zero quasiparticles, the 
second for two (at energies Ei and Ej ), and additional 
terms for four, six, etc. Z 0 similarly runs over odd oc¬ 
cupied states. Rewriting these sums as integrals over 
positive energies yields 

pOO 

F 0 — F e = —AiBT’lntanh / d E D(E) lncoth(E/2kBT), 

J o 

( 6 ) 


where D(E) is the density of states of the HQD, 

D(E)=p BGS (E) + ±p+(E) + lpo(E). (7) 
We take Pbgs{ e ) to be a standard BCS density of states, 


Pbcs(E) 


PaiVE 

y/E 2 - A (B) 2 


0(E - A) 


( 8 ) 


(6 is the step function), and po to be a pair of Lorentzian- 
broadened spinful levels symmetric about zero, 


Po( E ) 


7/27T 

(E-E±) 2 + ( 7 /2) 2 


7/2-7T 

(E + E±) 2 + ( 7 /2) 2 ' 


(9) 

Zeeman splitting of the bound state and pair-breaking by 
the external magnetic field are modeled with the equa¬ 
tions 


E t( B ) = A{ f } E 0 ±l 2 gp B B, ( 10 ) 

a(B > = a \M£) 2 - < n > 


where Eq is the zero-field state energy and A is the zero 
field superconducting gap. In the event that a bound 
state goes above the continuum, E+ > A (B), we no 
longer include the state in the free energy. Equation ([6]) 
was integrated numerically to obtain theory curves in 
Fig. (4). 

Equations (101 and © are reasonable provided the 
lower spin-split state remains at positive energy, E q > 0. 
For sufficiently large B c , the bound state will reach zero 
energy, resulting in topological superconductivity and 
Majorana modes, the subject of future work. 
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SUPPLEMENTARY INFORMATION: 

QUASIPARTICLE PARITY DYNAMICS IN A SUPERCONDUCTOR-SEMICONDUCTOR HYBRID 

QUANTUM DOT 


1. Effective gate voltage definition 

2. Bound on the single quasiparticle relaxation time 

3. Detailed interpretation of Coulomb diamonds 

4. Derivation of transport and thermal model 

5. Comparison of free energy approximations 

6 . Effect of the bound state on the free energy 


[TJ Effective gate voltage definition 

We define an effective gate voltage in software to tune the HQD. The physical gate voltages, Ur and Vl, are related 
to the effective gate voltage, Uq, by 


Ur — Vr, 0 + K ^g, 

U l = U L , 0 + Vl-K 2 U G , 

with k = 0.9997 and offset voltages Vr,q = —2.41 V, Vl,o = —3.96 V. 

These transformation rules ensure that Uj = (Vr — Vr,q) 2 + ( Vl — Vl,o) 2 , so that Vq can be interpreted as the 
distance from (Vr, 0 , Vl,o) in the Vr — Vl plane. 

All measurements are performed at backgate voltage Ubg = 2.39 V. 


[2] Bound on the single quasiparticle relaxation time 


The effect of quasiparticle relaxation is shown in Figs. Sla-e. Quasiparticle relaxation results in a disappearance 
of the negative differential conductance, in combination with the appearance of an extra conductance threshold. We 
quantify this observation by introducing the relative conductance ratio 


R = g' + ffNDC 
g' - ffNDC 


(SI) 


where <?ndc is the minimum of the negative differential conductance, and g' is the maximum of the extra conductance 
threshold that appears when r qp -A 0 (see Fig. Sli). The i?-value is a metric for the relative strength of the negative 
differential conductance. 

Figs. Slf-j show example conductance traces at constant bias and their associated f?-values. The traces show that 
R « — 1 corresponds to slow quasiparticle relaxation, and R ss 1 corresponds to fast quasiparticle relaxation. 

Fig. S2 shows the i?-value calculated as a function of single quasiparticle relaxation time, r qp . Also shown is a 
measured i?-value averaged over all negative differential conductance features in Fig. 1 of the main text. The measured 
R -value is consistent with r qp >0.1 /zs, giving the experimental bound on the single quasiparticle relaxation time. 


[3| Detailed interpretation of Coulomb diamonds 

Each conductance threshold in the Coulomb diamond plots can be interpreted with the help of the transport model, 
as shown in Fig. S3. For example, the highest bias at which NDC is observed occurs at the intersection of black and 
green lines, when vsd = (A + Eo)/e. 
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Fig. SI: Effect of quasiparticle relaxation, a Measured conductance g versus source-drain bias Vsd and gate Vg- b, 
Transport model of a, with r qp = oo. vg = uVg up to an offset, where a is the gate lever arm. Axis units are A/e = 180 gV. 
c-e, Model with r qp = 0.1 gs, r qp = 5 ns, and r qp = 0 respectively, f-j, Conductance versus gate at constant bias indicated in 
a. Relative conductance ratio, R = (</ + <?ndc)/(</ — <7ndc), for theory curves is labeled (see text). 



Fig. S2: Single quasiparticle relaxation bound. Relative conductance ratio, R = (g 1 + gtmc)/(g' — <?ndc), versus single 
quasiparticle relaxation time r qp . Dashed curve is theory derived as shown in Fig. SI. Data is the average over all charge 
transitions in Fig. 1, with vertical error the standard deviation of the mean, and horizontal error propagated from vertical. 


[4j Derivation of transport and thermal model 

This section gives a detailed derivation of the transport and thermal model used in the main text. 

To describe the electron transport through a metallic superconducting quantum dot we consider the following 
model: 


H = H hR + H u + H Tl 


(S2) 
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Fig. S3: Interpreting conductance thresholds, a, Calculated conductance g versus usd and vg- Eq = 0.1 A, all 
other model parameters same as main text. Dotted lines are dsd /2 = ±(ug + A/e) [black], dsd /2 = ±(ug + Ea/e) [blue], 
^sd /2 = ±(vg — Eo/e) [green], vsd /2 = ±(vg — A/e) [red], b, Eq = 0.32A, all model parameters same as main text, c, 
Eq = 0.5A, all other model parameters same as main text. Color scale shared across all plots. 


where the Hamiltonian 


Hlr = 53 (e 

OLVS 


Oil'S 


E a )< 


'cut's on's 


(S3) 


describes the normal metallic leads with cj. tvs being an electron creation operator in the lead a £ {L,R}, with an 
orbital quantum number v and spin s £ {t^ 4-}- The leads have chemical potentials given by /r Q = ±R/ 2 , where V 
denotes symmetrically applied bias. For the semiconductor-superconductor hybrid quantum dot, we use a simplified 
model consisting of a Bardeen-Cooper-Schrieffer (BCS)fSll continuum and an Andreev bound state for fixed number 
of particles [S2j : 


Hu 


^ "j A ^ ^ H n ^y nS p r y nS p 

s,p=e,h n 


■E 


TV) 


(S4) 


E c n = U(N - N g f , N = 4 s d 0s + ^ did, 

s n 


(S5) 


where d' ns creates an electron in the continuum with quantum number n (e.g. momentum of electrons on the dot), 
and dl denotes electron creation in a localized level, which gives rise to a subgap state in the BCS spectrum. The 
charging effects on the quantum dot are described by constant interaction model given by the term i?^-, where the 
charging energy is given by Eq = 217 and the number of electrons on the dot is controlled by a gate voltage which 
is parameterized by dimensionless number N g = V g /Ec■ The operator N gives the total number of electrons on the 
dot. With superconducting pairing, the dot Hamiltonian is diagonal in the basis of the quasiparticle operators 7 3 
which are given by [S3HS5| 


d ns = u nlns,e + SV nl- n s, h , 

(S 6 a) 

7 ns,e = U nd ns ~ SV n d_ ng S, 

(S 6 b) 

lns,h = u nSU ns - SV n dl nS , 

(S 6 c) 

^ns^e 7ns ,hi ^ns,h ^ 7ns,e) 

(S 6 d) 
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with denoting the Cooper pair creation operator and 7 ^ s e , h denoting the quasiparticle creation operator, which 
adds an electron/hole to the system. Here a state with quantum numbers — ns is the time-reversed partner of a state 
with quantum numbers ns. The quasiparticle excitation energy E n and the BCS coherence factors u n , v n are given 
in terms of superconducting gap A and electron dispersion on the dot e n as 

E ~ = En At = = <S7) 

Similarly, the subgap state operator is expressed as 

d 0s = u Ol 0 s,e + sv o7os,h > ( S8 ) 

with u o, Vq being model dependent coherence factors, which we set to uq = Vq = l/v^2 for simplicity. Lastly, the 
electrons in the leads and in the dot are coupled by the tunneling 

^ ’ ^S,aC Q , I/s dQ s T* ts,a d 0s^a^s T ^ . ( J'C,cx d a l , s d ns tQ a d ns C al/ ^j 

on's n 

= [*S,acL s K7 0s ,e + Sv 07og,h) + f S,a( U o7o s ,e + SV o7 0 s,h) C co,s (S9) 

Oil'S 

+ {fe, a cL s (M„ s ,e + SV nl-ns,h) + *C,a( U nT L,e + SV u7- nS ,h) C a ^}\ , 
n 

where tc, a gives the tunneling amplitude to the continuum and ts, a gives the tunneling amplitude to the subgap 
state. 


Thermodynamics of the even/odd effect 


We now present the free energy difference between the superconducting metallic island having even or odd number 
of electrons. The parity of the number of quasiparticles has to be equal to the parity of the number of electrons N 
on the island. The free energy difference SF between the odd and even occupation is expressed as 


«P = P 0 -P. = -il.(| 

in terms of the partition functions for different parities 

2 Z o/e = [](1 + e-^") T I^ 1 - e_/3B ")‘ 


(S10) 


(Sll) 


where /3 = 1 //cbT denotes the inverse temperature of the island. For a sufficiently large island the single particle 
spectrum can be described by the spectrum of a grounded superconductor where the single particle spectrum E n is 
given by Eq. 


Without a subgap state the free energy difference Eq. (|S10|) is expressed as 

i£lncotl/^ 1 


^-Fbcs = lntanh 


0E n 


lntanh J dE pbcs{E) lncoth > (S12) 


where pbcs(-E') is the BCS density of states for quasiparticles on the island given by 

Pbcs{e) = Ve^-a^ (S13) 

with pb = Pa\V denoting the normal state density of states, including spin, and pai is aluminum density of states per 
volume, and V is the volume of the island. For small temperatures /3A 1, the free energy difference (S12) can be 

approximated as 


SFbcs ~ —k&T lntanh 


r+oo 


dEpBcs{E)e dE 


= —UbT lntanh [A e ge /3A ] « A — fcBTln(A e ff), 


(S14) 
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where an effective number of quasiparticle states N e g is given by 

/»+0O 


r~ l-oo 

N eS = 2 dEp BCS {E)e-^ E - A ^ = 2p B >Ae /3A K 1 (/3A) « p D ^2^TA 

J A 

and K w (x) denotes the modified Bessel function of the second kind. 

With a subgap state the free energy difference Eq. (S10) acquires an additional term and one gets 

OO 


(S15) 


(5-Fabs = — fcsTlntanh 


U A 


dE pegs (E) In coth ( ^ ) + lncotli(/3£ , 0 /2) 


(S16) 


See also Eq. (S421 where the approximate expression for the first term in used. 

In the main text we discuss how the low-temperature data deviates from the above Andreev-bound-state model in 
terms of a life-time broadening of the subgap state. This is done by including a phenomenological broadening with 
width 7 into the subgap density of states, which then gives the free energy difference 


^'Aabs-lb = — /csTlntanh 


P+OO 


dE pscs {E) In coth ( ) + b ^2 


1 


T-± 1 1 

S=td 


r+°°dw 7 In coth 
j 2 tt (w - tE 0 ) 2 + (7/2) 2 


(S17) 


In the kinetic equation calculation presented below, the equilibrium distributions of quasiparticle in the continuum with 
an even or odd number of quasiparticles are needed. Since we will assume that the particles occupying the continuum 
are effectively equilibrated, we find the distribution functions by modifying the usual Fermi-Dirac distribution function 
as 

UE) = ^“ + 1> (S18) 

= e P(E-6F BC s ) + 1’ 

where P £ {e,o}, and P represents the opposite of P. 


fp(E) = 


eP E {Z P /Zp) + 1 


Number of quasiparticles 


Using the above results, we derive a simple expression for the number of quasiparticles in the absence of a bound 
state. At low temperature, when 5F B cs = A — fcBEln(A e ff), the distribution functions take the form 


fe 

fo 


N eS e~^ E+A \ 
1 C -/3(E- A)^ 


(519) 

(520) 


where N e g is given by Eq. (S15). 

The number of quasiparticles in each parity state can then be calculated using 

/> + OO 


N P = 2 


d E PBCs(E)fp(E). 


(S21) 


Substituting the above expression for f 0 gives N a = 1, as expected. Substituting f e gives the quasiparticle number 

/» + OO 


N e = 2 I dE p BC s(E)N eS e-^ E+A) 

J A 

r- l-oo 

= N eS e~ 2l3A / dE 2p BC s(E)e-^ E - A ^ 
J A 

= (iVeffe-^) 2 . 


(S22) 


Because of the large charging energy N e is the square of the bulk value N e ge 13A , indicating that quasiparticles must 
be created in pairs. 
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Incoherent transport 


We now calculate the current through the quantum dot by a set of master equations, with transition rates calculated 
by Fermi’s golden rule, which is valid in the weak tunneling regime or the so-called sequential tunneling regime. 

According to the Fermi’s golden rule the transition rate from the initial state i to the final state j caused by a 
perturbation V is given by |S8) : 


F/» = 2n\(f\V\i)\ 2 6{E f - Ei). 


(S23) 


In our case we are interested in the transitions caused by the tunneling Hamiltonian (S9| between the states "D) |LR) 
and |£>')|LR'): 


IW, L r'lr = 27t|(LR'| {D'\Ht |2?)|LR) \ 2 S(Ep> - E v + FIlr' - E LR ), (S24) 


where the many-body eigenstates of the lead Hamiltonian F/lr are denoted as |LR) and of the dot Hamiltonian FFd 
are denoted as \T>). Also Ep gives the energy of the state \V) and F? L r gives the energy of the state |LR). So we see 
that we need the following matrix element 

|(LR'|(P'|FF T |P)|LR)| 2 = Y |t a | 2 {| U „| 2 |(LR'|cL s |LR)| 2 |(2? / |7„ s , e |©)| 2 

otvsn 

+KI 2 | (LR' |cL s |LR) | 2 1 (V W_ n - Bih \V) | 2 (S25) 

+| Url | 2 |(LR'| Cal/s |LR)| 2 |(I?'| 7 t Sie |2?)| 2 

+K | 2 1 (LR' |c Q „ s |LR) | 2 1 (£> , |7_„ v jr>)| 2 } • 

Note that we have not written out the terms with the subgap state, but they have an analogous structure. 

In order to form the tunneling rates we need to specify what kind of states \D) on the metallic superconducting 
quantum dot we consider and what kind of (thermal) averaging procedure for these states we employ. We classify the 
states \T>) according to the number of electrons IV, the number of excitations Ns in the continuum, quantum number 
l labeling the configuration of the excitations in the state, and the occupancy of the subgap state x £ 2}, i.e., 

\V) = \N,N s ,l,x). (S26) 


The above states have the energy 


N s 

E N>Ns , l}X =E c N + E x + Y, E n■ (S27) 

nEZ 

If the number of electrons N is even/odd then the total number of excitations N to t,s also has to be even/ood. This 
means that the number of excitations in the continuum Ng depends on the subgap state occupancy x and the charge 
on the dot N, i.e., 


N s 


even for N even and x £ {0, 2}, 
odd for N even and x € {t, I}, 
odd for N odd and x £ {0, 2}, 
even for N odd and x £ {f, !}■ 


(S28) 


Now we want to thermally average over all quasiparticle states Ng of the continuum and their configurations l. 
The thermal averaging assumption is valid if the relaxation rate of quasiparticles is much faster than the rate of the 


tunneling events between the leads and the island. From the expression (S25) we see that we need to consider the 
thermal averages of the following expressions for the dot 


fs(E n ,N,x) = Y W n ,n s ,i,x\( N ’ N s,h x\'lL,p'yns,p\N, A' s ,Z,x)| 2 , 

l,N s 

f s (E n ,N,x) = Y W n> n s ,i,x\{ n , N s,l,x\'y n s lP 'yL, P \ N > N s ,l,x)\ 2 = 1 - f s {E n ,N,x ), 

1,N S 


(S29a) 


(S29b) 
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where p = e, h. Here Wn,n s ,1 denotes a thermal distribution for which we have 


Wn,N s ,1,x 


_ e — PE N ,N S ,l,: 

r. 


y wn,n s ,i,x- 

Ns,l 


By using Eq. (S18) and following the prescription (S28), we get the distributions 


{ fe(E), for N even and x = 0, 2, 
fo(E), for N even and x =t,l, 
fo(E), for N odd and x = 0, 2, 
f e (E), for N odd and x =t,4~ 


(S30) 


(S31) 


After thermally averaging over the source-drain lead states |LR), using grand-canonical ensemble, and the continuum 
states of the dot, we obtain the following tunneling rates from and to the continuum of the dot 


r+°° EdE 

r Q AT+ x <-AT ~ la / , = f¥i{E C N+x — E C n + x'E — XHa)fs,-x’(E, N, x ), 

Ns+x'^Ns J\A\ V E - l A l 

with fs,+ {E, N,x) = fs{E,N,x), /s ,-(E, N, x) = f s (E, N, x). 


(S32) 


Here = 2 x 27t x PBPa\tc,a\ 2 with p a denoting density of states of the normal leads. The continuum coupling 
7 a = 2 x 27r x PT>Pa\tc,a\ 2 is related to the normal-state conductance by gA\ = (n/2)(e 2 /h)"f a . Note that we have set 
the chemical potential of the metallic superconducting dot at zero go = 0 in order not to complicate the calculations, 
and also used that e n = 

Additionally, there are tunneling rates from and to the subgap state. When the starting state has no quasiparticles 
in the subgap state, i.e., |0), we get the following rates 


^n-u-n ~ r a Uo[l — /n(^at — E s — E^ r _ 1 — g a )], (S33a) 

Si — 0 

^n+u-n ~ r a ulh(E c N+1 + E s — E 1 ^ — p a ). (S33b) 

Si — 0 

For the state with single quasiparticle |s) we get 

Tat-k-at ~ r««Q[l — /n(^at + E s — E c N _ l - p a )\, (S34a) 

0 «-s 

TaT-K-AT ~ — /n(^AT — Eg — E C n _ 1 — g a )], (S34b) 

2 «-s 

^n+u-n ~ r Q iiQ/N(.Ejv + i + E s — Ex — g a ), (S34c) 

2 i—S 

^n+u^n ~ ra^o/N^Ar+i — E s — E% - g a ), (S34d) 

0 ^—s 

and for the state with two quasiparticles |2) we get 

Tat-k-at ~ r««Q[l — /n(^at + E s — E c N _ l - p a )], (S35a) 

si- 2 

^N+l<-N ~ ^aVof^i^N+l ~ Eg — E C N — p a ). (S35b) 

si-2 


Here T Q = 27rp a |ts,a| 2 , and s £ {ti4-} with s denoting the opposite of s. We include the relaxation from the continuum 
to the subgap state by introducing the following rates within the same charge state 


r N 0 i—N 0 = Tjv e <-Ar e = r relax . (S36) 

0-<—S 2 i —s 

Now we want to find the current through the superconducting metallic quantum dot. To do this we need to obtain 
the occupation probabilities Pn,x of the states described by a number of electrons N on the dot and the occupancy 
of the subgap state x. We write the following steady state Pauli master equation for the probabilities Pn,x 


d t Pn ’ x ~ 


£ 

N'x' 


r 


A7'<— A7-' N,x 


£ 

N'x' 


Ni^N 


Pn'x' = 0, 


(S37) 






with the condition 


E Pn,x = 1- 

N,x 


The rates entering in (S371 are given by 


"v _ 'pL _i_ pi? 

- N^i—N — 1 N'<-N "r" 1 N'<-N ’ 
■<—a: x'i—x x'-i—x 


and for x = x' we have 


pa _ pQ,2) 

1 N'*-N ~ 1 N'l—N 
x<—x 


N'^N • 
Ns — !•<— iVs JVs+l-<—iVs 


When the occupation probabilities Pn,x are obtained, the current through the quantum dot can be written as 


heq = (—e) £ ( r 


N,c 


■ N+U-N ~ ^N-li-N I r N,x- 
x'<—x x'<—x 


Pa 


(538) 

(539) 

(540) 

(541) 


These formulae are then solved numerically to produce the plots in Fig. SI and Fig. 2b in the main text. 


|5j Comparison of free energy approximations 


This section gives examples of the free energy difference, F a — F e , calculated under different approximations, 
considering the case without broadening 7 = 0 and without an applied field B = 0. Under these conditions the free 
energy difference is given by Eq. ( |S16| ). When (3 A > 1 the approximation lncoth(/3P/2) ss 2e~ l3E can be used for the 
first term. Applying the identity J A °° dE pBCs(E)e~ l3E = PaiVAK i(/3A) then gives 


F - F 

1 o 1 e 


—feeTlntanh 


2p^\V AKi(/3A) + lncoth 



(S42) 


where Ki{x) is a Bessel function of the second kind. In the very low temperature limit j3A 1, (3Eq > 1 the 
approximations Ki{f3A) « ^/ 7 r/( 2 / 3 A)e _/3A , lncotli(/3U 0 /2) « 2e~^ E °, and tanh(a:) « x can be used, giving 


F a -F e 


/cbPIh 


N eS e~ fiA + 2e~ 0Eo 


(S43) 


where N e g = pa\V y/2nk^T A. 

Equations (S42) and (S43) constitute two levels of accuracy at which Eq. (S16) can be evaluated. Figure S4 
compares the methods. Equation (S42) is an excellent approximation to Eq. (S16l over the experimentally relevant 
temperature range. Equation (S43) is poor approximation at intermediate temperatures. 


[b] Effect of the bound state on the free energy 

Figure S5 shows a comparison of the free energy difference, F 0 — F, with and without the subgap bound state. 
As the lowest energy unoccupied state, the bound state causes the free energy to saturate at F a — F e = E 0 at low 
temperature. It should be noted that the free energy difference with a subgap bound state was also shown in Fig. 5 
of Lafarge et al. fS7j . 


[51] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Phys. Rev. 108, 1175 (1957) 

[52] M. Tinkham, Phys. Rev. B 6, 1747 (1972) 

[53] B. D. Josephson, Phys. Lett. 1, 251 (1962) 

[54] J. Bardeen, Phys. Rev. Lett. 9, 147 (1962) 

[55] M. Tinkham, Introduction to Superconductivity (Dover Publications, 2004). 

[56] M. Tuominen, J. Hergenrother, T. Tighe, and M. Tinkham, Phys. Rev. Lett. 69, 1997 (1992) 

[57] P. Lafarge, P. Joyez, D. Esteve, C. Urbina, and M. Devoret, Phys. Rev. Lett. 70, 994 (1993) 

[58] H. Bruus and K. Flensberg, Many-Body Quantum Theory in Condensed Matter Physics (Oxford University Press, 2004). 


























9 



Fig. S4: Comparison of free energy approximations. Free energy difference F 0 — F e versus temperature 
different expressions for the free energy. All parameters same as main text (A = 180 /teV, Eq = 58 /teV, 7 = 0, B 
crosses are numerically exact values from Eq. (S16 1 , red line is Eq. (S421, blue line is Eq. (S43|. 


T for three 
= 0). Black 



Fig. S5: Effect of bound state on free energy. Free energy difference F a — F e versus temperature T with and without the 
semiconductor bound state. All parameters the same as main text (A = 180 /mV, Eq = 58 /reV, 7 = 0, B = 0). Blue trace 
includes the BCS density of states only, red trace includes the BCS density of states and the discrete state. 









